###########################################################################
## The code below replicates the figures and tables in:
## Sharan Grewal, "Norm Diffusion through US Military Training in Tunisia,"
## Security Studies, forthcoming. 
## The code requires two datasets: "survey1.csv" and "survey2.csv"
###########################################################################

## Load packages
library(stargazer)
library(ggplot2)
library(scales)
library(effects)

## Load data
data <- read.csv("survey1.csv")
data2 <- read.csv("survey2.csv")

## Set reference group for foreign training
data$train <- factor(data$train, levels=c("Only France","US","Other"))
data2$train <- factor(data2$train, levels=c("Only France","US","Other West","Other","No For Train"))



##########################
## Figure 2: Histograms ##
##########################

## Survey 1: Right to Vote
h = hist(data$vote, breaks=0:5)
h$density = round(h$counts/sum(h$counts)*100, 1)
plot(h, freq=FALSE, col="grey", labels=TRUE, ylim=c(0, 100),
     xaxt='n', xlab=NA, ylab="Percent",
     main="Military should have right to vote (N=70)")
axis(1, 0.5:4.5, c("Strongly\nDisagree", "\nDisagree", "\nNeutral", "\nAgree", "Strongly\nAgree"), tick=FALSE)


## Survey 1: Retired Officer as MOD
h = hist(data$retiredMOD, breaks=0:5)
h$density = round(h$counts/sum(h$counts)*100, 1)
plot(h, freq=FALSE, col="grey", labels=TRUE, ylim=c(0, 100),
     xaxt='n', xlab=NA, ylab="Percent",
     main="Defense minister should be retired officer (N=70)")
axis(1, 0.5:4.5, c("Strongly\nDisagree", "\nDisagree", "\nNeutral", "\nAgree", "Strongly\nAgree"), tick=FALSE)


## Survey 2: Voting
h = hist(data2$voting, breaks=0:5)
h$density = round(h$counts/sum(h$counts)*100, 1)
plot(h, freq=FALSE, col="grey", labels=TRUE, ylim=c(0, 100),
     xaxt='n', xlab=NA, ylab="Percent",
     main="Voting in Elections is... (N=159)")
axis(1, 0.5:4.5, c("Very\nInappropriate", "\nInappropriate", "\nNeutral", 
      "\nAppropriate", "Very\nAppropriate"), tick=FALSE, cex.axis=0.9)


## Survey 2: Retired Officer as President
h = hist(data2$retiredpres, breaks=0:5)
h$density = round(h$counts/sum(h$counts)*100, 1)
plot(h, freq=FALSE, col="grey", labels=TRUE, ylim=c(0, 100),
     xaxt='n', xlab=NA, ylab="Percent",
     main="President Should be a Retired Officer (N=190)")
axis(1, 0.5:4.5, c("Strongly\nDisagree", "\nDisagree", "\nNeutral", 
                 "\nAgree", "Strongly\nAgree"), tick=FALSE)



#############
## Table 1 ##
#############

one <- lm(retiredMOD~train+sex+army+interior+MOD+essebsi+democracy+ryear+factor(rank), data=data)

two <- lm(vote~train+sex+army+interior+MOD+essebsi+democracy+ryear+factor(rank), data=data)

three <- lm(retiredpres~train+sex+army+interior+MOD+essebsi+democracy+active+factor(rank), data=data2)

four <- lm(voting~train+sex+army+interior+MOD+essebsi+democracy+active+factor(rank), data=data2)

stargazer(one, two, three, four, 
      dep.var.labels=c("Retired MOD","Right to Vote","Retired Pres","Voting"),
      covariate.labels=c("U.S. Training", "Other Western Training","Non-Western Training","No Foreign Training","Female","Army", "Interior Region", "MOD Horchani/Zbidi", "Pres. Essebsi","Democracy", "Retirement Year", "Major","Colonel","Colonel-Major","Active-Duty"))




##############
## Figure 3 ##
##############

## Survey 1: Right to Vote
data.summary <- predict(two, newdata=data.frame(train=c("Only France","US"), sex=rep(mean(data$sex),1), ryear=rep(mean(data$ryear),2), army=rep(mean(data$army),2), interior=rep(mean(data$interior),2), MOD=rep(mean(data$MOD),2), essebsi=rep(mean(data$essebsi),2), democracy=rep(mean(data$democracy),2), rank=rep(4,2)), interval="confidence", level=0.9)
data.summary <- data.frame(treatment=c("French Training (N=32)", "US Training (N=39)"), data.summary)

ggplot(data.summary, 
       aes(x=treatment, y=fit)) +  
  geom_bar(position=position_dodge(width=0.9), stat="identity", fill=c("darkgrey", "darkblue")) + 
  geom_errorbar(aes(ymin=lwr, ymax=upr, width=0.1),
                position=position_dodge(width=0.9)) +
  ggtitle("Right to Vote") +
  theme_bw() + 
  theme(panel.grid.major = element_blank()) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_y_continuous(limits=c(0.5, 5),breaks=seq(1, 5, by=1),oob=rescale_none) +
  xlab("") +
  ylab("Level of Agreement (1-5)") +
  geom_text(data=data.summary, 
            aes(x=treatment, y=fit,
                label=round(fit,1)),
            vjust=-2.5,size=7,position=position_dodge(0.9)) +
  theme(text = element_text(size=18)) +
  geom_text(aes(x=1.5, y=4, label="p=0.043"), size=7)


## Survey 1: Retired Officer as MOD
data.summary <- predict(one, newdata=data.frame(train=c("Only France","US"), sex=rep(mean(data$sex),1), ryear=rep(mean(data$ryear),2), army=rep(mean(data$army),2), interior=rep(mean(data$interior),2), MOD=rep(mean(data$MOD),2), essebsi=rep(mean(data$essebsi),2), democracy=rep(mean(data$democracy),2), rank=rep(4,2)), interval="confidence", level=0.9)
data.summary <- data.frame(treatment=c("French Training (N=32)", "US Training (N=39)"), data.summary)

ggplot(data.summary, 
       aes(x=treatment, y=fit)) +  
  geom_bar(position=position_dodge(width=0.9), stat="identity", fill=c("darkgrey", "darkblue")) + 
  geom_errorbar(aes(ymin=lwr, ymax=upr, width=0.1),
                position=position_dodge(width=0.9)) +
  ggtitle("Retired Officer as MOD") +
  theme_bw() + 
  theme(panel.grid.major = element_blank()) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_y_continuous(limits=c(0.5, 5),breaks=seq(1, 5, by=1),oob=rescale_none) +
  xlab("") +
  ylab("Level of Agreement (1-5)") +
  geom_text(data=data.summary, 
            aes(x=treatment, y=fit,
            label=round(fit,1)),
             vjust=-2.5,size=7,position=position_dodge(0.9)) +
  theme(text = element_text(size=18)) +
  geom_text(aes(x=1.5, y=4, label="p=0.006"), size=7)



## Survey 2: Voting
data.summary <- predict(four, newdata=data.frame(train=c("No For Train","Only France","US","Other West"), sex=rep(mean(data2$sex==1),4), active=rep(mean(data2$active, na.rm=T),4), army=rep(mean(data2$army==1),4), interior=rep(mean(data2$interior),4), MOD=rep(mean(data2$MOD),4), essebsi=rep(mean(as.numeric(data2$essebsi),na.rm=T),4), democracy=rep(mean(as.numeric(data2$democracy), na.rm=T),4), rank=rep("Officer",4)), interval="confidence", level=0.9)
data.summary <- data.frame(treatment=c("No Foreign\nTraining (N=187)", "French Training\n(N=19)","US Training\n(N=28)", "Other West\n(N=13)"), data.summary)
data.summary$treatment <- factor(data.summary$treatment, levels=c("No Foreign\nTraining (N=187)", "French Training\n(N=19)","US Training\n(N=28)", "Other West\n(N=13)"))

ggplot(data.summary, 
       aes(x=treatment, y=fit)) +  
  geom_bar(position=position_dodge(width=0.9), stat="identity", fill=c("lightgrey", "darkgrey","darkblue","grey")) + 
  geom_errorbar(aes(ymin=lwr, ymax=upr, width=0.1),
                position=position_dodge(width=0.9)) +
  ggtitle("Voting is appropriate") +
  theme_bw() + 
  theme(panel.grid.major = element_blank()) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_y_continuous(limits=c(0.5, 5.1),breaks=seq(1, 5, by=1),oob=rescale_none) +
  xlab("") +
  ylab("Level of Agreement (1-5)") +
  geom_text(data=data.summary, 
            aes(x=treatment, y=fit,
                label=round(fit,1)),
            vjust=c(-2.1,-3.1,-2.9,-3.4),size=7,position=position_dodge(0.9)) +
  theme(text = element_text(size=18)) +
  geom_text(aes(x=2.5, y=5, label="p=0.046"), size=7)



## Survey 2: Retired Officer as President
data.summary <- predict(three, newdata=data.frame(train=c("No For Train","Only France","US","Other West"), sex=rep(mean(data2$sex==1),4), active=rep(mean(data2$active, na.rm=T),4), army=rep(mean(data2$army==1),4), interior=rep(mean(data2$interior),4), MOD=rep(mean(data2$MOD),4), essebsi=rep(mean(as.numeric(data2$essebsi),na.rm=T),4), democracy=rep(mean(as.numeric(data2$democracy), na.rm=T),4), rank=rep("Officer",4)), interval="confidence", level=0.9)
data.summary <- data.frame(treatment=c("No Foreign\nTraining (N=187)", "French Training\n(N=19)","US Training\n(N=28)", "Other West\n(N=13)"), data.summary)
data.summary$treatment <- factor(data.summary$treatment, levels=c("No Foreign\nTraining (N=187)", "French Training\n(N=19)","US Training\n(N=28)", "Other West\n(N=13)"))

ggplot(data.summary, 
       aes(x=treatment, y=fit)) +  
  geom_bar(position=position_dodge(width=0.9), stat="identity", fill=c("lightgrey", "darkgrey","darkblue","grey")) + 
  geom_errorbar(aes(ymin=lwr, ymax=upr, width=0.1),
                position=position_dodge(width=0.9)) +
  ggtitle("Retired Officer as President") +
  theme_bw() + 
  theme(panel.grid.major = element_blank()) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_y_continuous(limits=c(0.5, 5),breaks=seq(1, 5, by=1),oob=rescale_none) +
  xlab("") +
  ylab("Level of Agreement (1-5)") +
  geom_text(data=data.summary, 
            aes(x=treatment, y=fit,
            label=round(fit,1)),
            vjust=c(-2,-2.6,-2.2,-3),size=7,position=position_dodge(0.9)) +
  theme(text = element_text(size=18)) +
  geom_text(aes(x=2.5, y=4.9, label="p=0.025"), size=7)




#############
## Table 2 ##
#############

one <- lm(democracy~train+sex+army+interior+MOD+essebsi+ryear+factor(rank), data=data)
two <- lm(chief~train+sex+army+interior+MOD+essebsi+democracy+ryear+factor(rank), data=data)
three <- lm(parli~train+sex+army+interior+MOD+essebsi+democracy+ryear+factor(rank), data=data)

four <- lm(democracy~train+sex+army+interior+MOD+essebsi+active+factor(rank), data=data2)
five <- lm(chief~train+sex+army+interior+MOD+essebsi+democracy+active+factor(rank), data=data2)
six <- lm(parli~train+sex+army+interior+MOD+essebsi+democracy+active+factor(rank), data=data2)

stargazer(one, two, three, four, five, six, 
      dep.var.labels=c("Democracy","Chief of Staff","Parli Advisor","Democracy","Chief of Staff","Parli Advisor"),
      covariate.labels=c("U.S. Training", "Other Western Training","Non-Western Training","No Foreign Training","Female","Army", "Interior Region", "MOD Horchani/Zbidi", "Pres. Essebsi","Democracy", "Retirement Year", "Major","Colonel","Colonel-Major","Active-Duty"))



## For questions or comments, please contact ssgrewal@wm.edu


